Numerical Study of a Superconducting Glass Model 
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An XY model with random phase shifts as a model for a superconducting glass is studied in two 
and three dimensions by a zero temperature domain wall renormalization group which allows one 
to follow the flows of both the coupling constant and the disorder strength with increasing length 
scale. Weak disorder is found to be marginal in two and probably irrelevant in three dimensions. 
For strong disorder the flow is towards a non-superconducting gauge glass fixed point in 2d and 
a superconducting glass in 3d. Our results are in agreement with recent analytic theory and are 
inconsistent with earlier predictions of a re-entrant transition to a disordered phase at very low 
temperature and with the loss of superconductivity for any finite amount of disorder. 
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The classical XY model with random phase shifts de- 
scribed by the Hamiltonian ||l| 



h = -J2j^^ 



9 J Aij) 



<ij> 



has been the subject of much interest over the past 
decade. Here, 9i is the phase of the order parameter at 
the i*'' site of a square lattice in 2d or a simple cubic lat- 
tice in 3d and the sum is over all nearest neighbor bonds. 
The coupling constants Jy are uniform Jy = J > 
and the quenched random phase shifts Aij are uncorre- 
lated from bond to bond and uniformly distributed over 
the range —air < Aij < air with < a < 1 so that 
< Aij >— and < \Aij\ >= a7r/2 where <> means an 
average over disorder. In d = 2, the model of eq.(P is a 
description of an array of Josephson junctions in a mag- 
netic field perpendicular to the plane of the array when 
9i represents the phase of the superconducting order pa- 
rameter of the «*'' grain and Aij = (27r/'I>o) Jr' A ■ dl with 
A the vector potential of the external magnetic field and 
$0 = hc/2e is the quantum of flux. The Aij become in- 
dependent quenched random variables when the average 
flux through an elementary plaquette is an integer mul- 
tiple of $0 but the superconducting grains are randomly 
displaced from their ideal lattice positions |g]. It is also 
a model for an XY magnet with random Dzyaloshinskii- 
Moriya interactions . 

Whatever the physical origin of the quenched random 
phase shifts, the system described by eq.(|l|) has been a 
theoretical challenge for a decade. Early work concluded 
that weak disorder (a <^ 1) does not destroy the super- 
conducting phase at intermediate temperature T but at 
low T there is a re-entrant transition to a normal phase 
[^,y|. However, this was not confirmed either numeri- 
cally nor experimentally B . Some later theoretical work 
Q suggested that the bound vortex (KT) phase [g| is de- 
stroyed at any finite T by arbitrarily small disorder and 
that the experimentally observed KT phase is a finite size 
effect. The more recent theoretical work [R-|l^l based on 
the ideas of Cha and Fertig |11|, on the other hand, argue 



for a more conventional phase diagram in which there is 
a superconducting phase for T < Tc{a) where Tc{a) > 
for a < ac- The case of maximum disorder (a — 1) when 



(1) the Aij are uniformly distributed beween — tt and 



is called the gauge glass [|12|,^ which has been studied 
numerically at T = by a domain wall renormalization 
group (DWRG) [|§^ in both d = 2 and d = 3 and 
the conclusions from these studies are that the stiffness 
to distortions of the phase vanishes in d = 2 so that the 
glass is not superconducting and that in d = 3 the gauge 
glass is probably superconducting, but the evidence is 
not conclusive. In view of the three conflicting scenarios 
in two dimensions, (i) re-entrant transition 00], (ii) de- 
struction of superconductivity for any finite disorder H] 
and (iii) superconductivity for T < Tc{a) with Tc{a) > 
for < a < Qfc 0,^,|l^, this system is an ideal candidate 
for study by a numerical DWRG at T = as this will 
distinguish scenario (iii) from the others as only scenario 
(iii) predicts a finite stiffness at T = for < a < ac. 

The standard DWRG @ at T = consists of com- 
puting the lowest energies of a set of systems of several 
linear sizes L with periodic and antiperiodic boundary 
conditions (BC) in one direction with some fixed BC in 
the other d — 1 directions. The difference AE{L) =< 
\Eap{L) — Ep{L)\ > is the domain wall energy and 
AE{L)/2 is interpreted as an effective coupling constant 
J{L) at length scale L which one expects to scale as 
J{L) ~ L^ at large L. The stiffness exponent is a cru- 
cial quantity as its value will distinguish between an or- 
dered superconducting phase at small but finite T {9 > 0) 
and a disordered phase {9 < 0). If 6* < 0, then the en- 
ergy of an excitation of size L is AE{L) which vanishes 
as L — > oo and the probability of such a phase unwinding 
excitation P{AE{L)) ~ exp{—AE{L)/kT) is large at any 
T > so the stiffness to twists in the phase will vanish. If 
9 > 0, the converse is true and the stiffness will be finite 
for T < Tc and the glass will be superconducting. This 
technique has been applied to random systems such as 
spin and gauge glasses at T = but the the value of the 
exponent 9 obtained by this method is not very reliable 
because it is not clear if the asymptotic scaling regime 



is reached for the small sizes L it is usually possible to 
simulate. 

For the problem of interest with a variable disorder 
strength, this version of the DWRG, which considers only 
the scaling of the effective coupling J{L) needs consider- 
able modification as we also want to know how the dis- 
order strength scales with L. For a single junction, the 
Hamiltonian is H = —Jcos{(j) — A) where (p is the phase 
difference across the junction and the usual comparision 
of the energies with periodic and antiperiodic BC gives 
AE{1) = 2JcosA which does not separate the disorder 
strength from the coupling constant. At length scale L, 
the interaction is Vl{<J) — A{L)) where A{L) is the phase 
shift at scale L and Vl {<p) is a 27r-periodic function with 
a minimum when its argument is zero. Thus, if we im- 
pose BC with phase shifts A^ across the boundaries in 
the d directions fi — 1,2, ...d, minimizing the energy with 
respect to the phases 6i will give the ground state energy 
of a system of linear size L as a function of A^ which is 
2-K periodic in each of the d directions 



El{^^) = El{^^ + 2tt) 



(2) 



with a minimum at some Aj^ which depends on the pre- 
cise realization of disorder. The key observation is that 
A° is exactly the phase shift A^(L) which minimizes the 
energy at scale L. A measure of the strength of disorder 
at this scale is 



\A{L)\ =< |A"| > 



(3) 



with |A(1)| — aTr/2. The coupling constant J{L) at scale 
L is found by first finding El{A^J, changing A" by n in 
one of the d directions and then finding the energy mini- 
mum El{A° + Tr) with these BC. As discussed above, the 
coupling constant J{L) at scale L is 



JiL)=<iELiA°+7r)-ELiA'))> 



(4) 



and measuring J{L) and |j4(L)| for several sizes L gives 
renormalization group flows for both the coupling con- 
stant and disorder strength. Of interest are the stable 
fixed point values J* = J{L = oo) and A* = \A{L = oo)| 
as these determine the nature of the phases. There are 
several possibilities of which the simplest are [J* — oo, 
A* = 0], [J* = oo. A* = 7r/2], [J* = 0, A* = tt/2] 
corresponding respectively to a superconducting state 
with long range order, a supeconducting glass and a non- 
superconducting glass. There are other possibilities such 
as a state with quasi long range order corresponding to 
a flow to a fixed line with finite J* and A* whose val- 
ues depend on the initial values of coupling and disorder. 
This is the scenario in d = 2 predicted by recent analytic 
work|,0. 

Of course, since we do not know how to find the exact 
ground state of a disordered system of arbitrary linear 
size L, the best we can do is to numerically estimate 
J{L) and |A(I/)| for a set of samples of different sizes 
L up to some maximum and extrapolate to large L. We 



use simulated annealing ||T^ to estimate the ground state 
energies which is considerably more efficient than simple 
repeated quenches to T = |1^. Also, we imposed peri- 
odic Ap = BC in d ~ 1 directions and twisted A y^ 
BC in the remaining direction and minimized the energy 
with respect to the the phases 9i and to the twist A to 
find A*'. To obtain the domain wall energy AEl the 
twist is changed to A*^ + n and kept fixed while the en- 
ergy is minimized with respect to the 9i only. According 
to our earlier discussion, the energy should be minimized 
with respect to global phase shifts in all d directions and 
the domain wall energy AE^ obtained by increasing the 
phase shift by n in one direction. To within the errors of 
our simulations, AE^ is independent of the choice of BC 
in the d — 1 transverse directions so, for simplicity, we 
imposed periodic or A^t = BC in these directions. As a 
consistency check |Q , we simulated two identical copies 
of each system with different random number sequences 
to obtain two estimates Ei , E2 of the ground state energy. 
In the event that the simulation finds the exact ground 
state, then 6E — Ei — E2 =0, which often ocurrs for 
our small L values. If the simulation does not reach the 
exact minima, < (SE)"^ > is a measure of the error. To 
minimize the errors caused by failure to reach the true 
energy minimum, we adjust the annealing schedule and 
the number of annealing attempts until SE/E < iV~^/^ 
where N — 10'^ in 2d and 10'* in 3d is the number of re- 
alizations of disorder. This consistency check makes the 
error due to not reaching the true ground state no worse 
than the statistical error in the averaging over disorder. 
To our knowledge, there is no analogue for this disor- 
dered XY system of the "branch and cut" algorithms 
to find exact ground states of Ising spin glasses |Q in 
fairly large systems in a reasonable amount of CPU time. 
For repeated simulated annealings of N different samples 
the CPU time becomes prohibitive for i > 8 in 2d and 
L > 4 in 3d. We therefore chose sizes L = 2,4, 8 in 2d 
and L = 2,3,4 in 3d and the results are summarized in 
Fig.(l) for 2d and in Fig. (2) for 3d. 

In 2d, for small disorder a < ac ~ 0.37, J{L) increases 
more slowly than a power of L and seems to flow to 
a finite disorder dependent value J* (a) and the disor- 
der strength |A(L)1 does not change with L, at least for 
our small sizes, both of which are completely consistent 
with analytic RG calculations [pUlO|]. At larger disorder 
strength a > ac, \A{L)\ increases and J{L) decreases as 
L increases, in agreement with the analytic theory p,p^. 
In this range of disorder strengths, the system is proba- 
bly flowing to the non-superconducting glass fixed point 
at J* — 0,A* — n/2. When the disorder is maximal 
(a = 1), |A(L)| remains fixed at tt/2 and J{L) ~ L^ 
with the stiffness exponent 9 w —1/2, in agreement with 
other simulations of the gauge glass in 2d |jl^,|l^ . The 
flows shown in Fig.(l) may be regarded as RG flows in a 
higher dimensional parameter space projected on to the 
(J(i), \A{L)\) plane, so the crossing of the two trajec- 
tories for a = 0.45 and a — 0.5 does not violate the 
non-crossing rule. For the 2d system, the RG flows are 



in the three parameter space of J, |A| and the vortex fu- 
gacity y ^j^^^. From the data of Fig.(l) and assuming 
that the trends for small L continue when L is large, 
one would conclude that the most likely scenario for the 
disordered system in 2d is, for weak disorder a < a^^ 
J{L) -^ J*{a) and \A{L)\ = |A(1)|, so that the cou- 
pling constant scales to a finite but disorder dependent 
value while the disorder is marginal. For larger disor- 
der a > ttc, J{L) -> J* = and \A{L)\ -^ tt/2 which 
is a non-superconducting disordered state. Our results 
are consistent in all respects with recent analytical the- 
ory pjl0| and inconsistent with the re-entrant |^,|[ and 
complete destruction of superconductivity p] scenarios. 
However, they are consistent with a re-entrant scenario 
in which the ordered phase extends to T = for a range 
of a pi]. The flows for 2d of Fig. (1) are consistent with a 
discontinuous jump in J* at ac as predicted analytically 



ponent 6 = 0. Although this is consistent with other 
simulations on the 3d gauge glass @,|6[ , our use of the 
phase representation of eq.(|lj) in the simulations together 
with the very small sizes may introduce large corrections 
to scaling. Nevertheless, the data is consistent with a 
superconducting glass phase which survives at finite T. 
These considerations suggest that the phase diagram for 
the model in 3d is similar to that of the corresponding 
infinite range model 123]. 
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FIG. 1. DWRG flows for the 2d superconducting glass 
model in the coupling constant-disorder strength (J, \A\) 
plane. The initial value of J is J(l) = 1 and the initial disor- 
der strength a is indicated near corresponding symbols. 

The results of the simulations of the random gauge 
model in 3d are shown in Fig. (2). The system sizes are 
very small {L = 2, 3, 4) because of the CPU time needed 
to get close to the minimum energies so that irrelevant 
variables are giving large corrections to scaling. Never- 
theless, some qualitative features are apparent, assum- 
ing that the small L trends continue. For small disorder 
a < ac ^ 0.55, J{L) ~ L'^~'^ as expected and the dis- 
order strength seems to decrease. It is impossible to say 
if |^(i)| ^ as one expects but the data is consistent 
with this. One is tempted to conclude that, in this regime 
of weak disorder, the DWRG flows are to a stable fixed 
point at J* = oo, A* = corresponding to a true su- 
perconducting phase. For larger disorder ac < a < n/2, 
the disorder increases with L and seems to flow to its 
maximum value of tt/2. The coupling J(L) seems to 
flow to a flnite value which corresponds to a stiffness ex- 
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FIG. 2. DWRG flows for the model in 3d. 

Our conclusions from our new T = DWRG which 
follows the flows in two parameter space are in 2d, the 
recent analytic theory which predicts a quasi long range 
ordered state for T < Tc{a) is the correct scenario and 
earlier suggestions of a re-entrant transition to a disor- 
dered phase at low T or no superconductivity at any fi- 
nite disorder are ruled out. In 3d, weak disorder has little 
or no effect on the supeconducting phase and there is a 
critical disorder strength parametrized by a = ac, above 
which the system is a superconducting glass at low T. 
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